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Turing instabilities for a two species reaction-diffusion systems is studied under anisotropic dif¬ 
fusion. More specifically, the diffusion constants which characterize the ability of the species to 
relocate in space are direction sensitive. Under this working hypothesis, the conditions for the onset 
of the instability are mathematically derived and numerically validated. Patterns which closely 
resemble those obtained in the classical context of isotropic diffusion, develop when the usual Tur¬ 
ing condition is violated, along one of the two accessible directions of migration. Remarkably, the 
instability can also set in when the activator diffuses faster than the inhibitor, along the direction 
for which the usual Turing conditions are not matched. 
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I. INTRODUCTION 

Spatio-temporal patterns are widespread in nature: beautiful spots and stripes appear on the coat of animals [9[, 
patterns of cracking emerge on the fracture surface of materials , reacting chemicals give rise to complex and dy¬ 
namical structures as in the celebrated Belousov-Zhabotinsky reaction spatial games in social sciences yield 

self-organized regular motifs sisiniiii- A common feature which is shared by the above mentioned applications 
is the spontaneous formation of complex structures, which result from the non trivial interplay between noise and 
deterministic dynamics. Elucidating the key mechanisms that seed the process of pattern formation is therefore an 
important topic of investigation of cross disciplinary impact. 

One of such mechanisms was identified and thoroughly discussed in a pioneering work of A. Turing [T5[: homo¬ 
geneous equilibrium solutions of a multi-species reaction-diffusion system can be destabilized upon injection of a 
small inhomogeneous perturbation. This latter undergoes an exponential amplification, in the linear regime of the 
evolution. Then, non linearities come into play and the system eventually reaches a patchy, spatially inhomogeneous, 
equilibrium. Traveling waves and spiraling patterns can be also generated following a Turing-like, symmetry breaking 
instability. 

In the classical setting, two mutually interacting species are considered: these are the so called activator and in¬ 
hibitor. If the diffusion is isotropic, or in other words not affected by the specific direction of displacement, the 
inhibitor species should diffuse faster than the activator, for Turing patterns to develop. Systems of three simulta¬ 
neously diffusing species [TsI have also been considered in the literature and shown to display a richer zoology of 
possible instabilities and pattern. In this generalized context, self-organized motifs can also develop, if one species is 
solely allowed to diffuse in the embedding medium Q. Beyond the deterministic scenario, stochastic Turing patterns 
have been also reported for reaction diffusion systems defined on a regular lattice or complex networks [l|, H, 0, 0 • 

Starting from these premises, and with reference to the paradigmatic scenario where just two species are made to 
interact, we shall here revisit the conditions that yield the Turing instability, under the assumption of anisotropic 
diffusion. More concretely, we shall derive sufficient conditions for the emergence of Turing like patterns in a rectan¬ 
gular, continuous, domain subject to periodic boundary conditions, assuming generic non linear reaction terms and 
imposing anisotropic, i.e. direction sensitive, diffusion coefficients. As we will demonstrate in the following, patterns 
do exist also if the condition for the onset of the Turing instability is uniquely satisfied along one direction. These 
latter patterns resemble quite closely those that are found under the standard assumption of isotropic diffusion, the 
non linearity being responsible for the mixing of cross modes. In addition, patterns can also flourish when the activator 
diffuses faster than the inhibitor, along one specific direction. In this case, the system organizes along the direction 


2 


orthogonal to the latter, hence displaying regular, just locally distorted, stripes. 

The paper is organized as follows. In Section [IT] we will present the reference framework and then, in Section Hill 
derive the mathematical conditions for the generalized anisotropic instability. Section [T^ is devoted to reporting some 
numerical tests to validate the theoretical analysis. Finally, we shall sum up and conclude. 


II. ANISOTROPIC DIFFUSION OF REACTIVE SPECIES ON CONTINUUM DOMAINS 


Let us consider two interacting species and denote by u and v their respective concentrations. The species can 
freely diffuse inside a rectangular domain, R = [0,1/^^] x [0,1/^] C M+ x M+, as specified by their respective diffusion 

coefficients. We shall in particular assume that the diffusion coefficients are anisotropic, meaning that they depend 

fee) 

on the specific direction of migration. More precisely, Du ^ > 0 denotes the diffusion coefficient for species u along 
direction x, while Du^ > 0 refers to the orthogonal direction y. Similar considerations respectively apply to > 0 

and > 0. The mutual evolution of species u and v is thus governed by the reaction diffusion equations: 


= /(w, v) + + D^y'^dlu 

= g{u, v) + + D^v'^dyV 


V(x, y) ^ R and Vt > 0 . 


( 1 ) 


where /(*,•) and are non linear functions of the concentration amounts. The above equations should be com¬ 

plemented by the initial conditions: 


u{x,y,0) = uo{x,y) and v{x,y,0) = vo{x,y) \/{x,y)eR, 


( 2 ) 


for some regular functions uq and vq, and suitable boundary conditions. In the following we shall adopt the Dirichlet 
periodic boundary conditions, namely 


u{x,0R) = u{x, LyR) VxG[0,La,] and 
u{Q,y,t) = u{Lx,y,t) My e[Q,Lv\ and 

and similarly for v. 

Let us assume the system m admits a stable, spatially homogeneous, 
translates in: 

such that: tr( J) = ^ < 0 and det( J) = fuQv - fvQu > 0 , (4) 

where J stands for the Jacobian matrix of system O: 

J=(fu fv\ ( 5 ) 

\9u 9v J 

where fu denotes the derivative of f{u,v) with respect to u, and similarly for fv^Qu^gv Here, and throughout the 
remaining part of the paper, we evaluate the partial derivatives at the equilibrium point {u^v). Without losing 
generality, we will also assume /« > 0 and < 0: is thus the activator species, while v refers to the population of 
inhibitors. 

The celebrated Turing patterns originate from a symmetry breaking instability of the homogeneous equilibrium 
solution. The introduction of an inhomogeneous perturbation around {u, v) activates the diffusion terms and, under 
specific conditions, makes the system to drift away from the deputed homogeneous equilibrium, towards a patchy, 
non homogeneous, asymptotically stable, solution. Mathematical conditions for the Turing instability to set in can be 
readily derived by first linearizing equations o and then Fourier transforming, both in time and space, the obtained 
linear system. This yields the so called dispersion relation, an equation for the growth rate Xk associated to Fourier 
mode k = {kx^ky). By carrying out this straightforward calculation, which is for instance detailed in Q, it can be 
eventually proven that Xk satisfies the following quadratic equations: 

+ d{kx^ ky)Xk + h(kx^ ky) = 0 , ( 6 ) 

where: 



Vt > 0 
Vt > 0 


( 3 ) 


solution u = u and v = v. This request 


dik,, ky) = -tr(J) + + Z)(-)) + kl{D^y) + ) 

ky) = det(J) - klifuD^:^ + g.Di^^) - + g.D^y^) 


( 7 ) 

(8) 
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Turing patterns materialize if the real part of takes positive values over finite window in k, which in turn amounts 
to require the presence of unstable non zero Fourier modes. We remark however that d(kx^ky) in Eq. ([6]) is always 

positive, since, by assumption, tr(J) < 0 and, in addition, > 0. Then, as a natural consequence, the Turing 

symmetry breaking instability can take place only if a compact domain exists in (kx^ky) such that h{kx,ky) < 0. As 
already mentioned, in the classical limit of isotropic diffusion, Du = D^f^ = and Dy = D^\ the Turing 

instability can take place only if the inhibitors diffuse faster than the activators, i.e. Dy > VcDu where Tc, the critical 
ratio of diffusivities, is a positive coefficient larger than 1. In the following we will show that this stringent assumption 
can be partially relaxed in the generalized setting where the diffusion constants are made to depend on the direction 
of propagation. 


III. TURING INSTABILITY IN PRESENCE OF ANISOTROPIC DIFFUSION 


The function h(kx^ky) is a multivariate polynomial of the variables k‘^ and ky. It is straightforward to check that 
it is positive at the origin and for large k‘^ and ky. We are here interested in determining when h{kx,ky) can change 
sign as function of k‘^ and ky^ so signaling the onset of the instability. To this end, we first consider restrictions of 
h{kx, ky) on kx = 0, and then on ky = 0. 

Focusing on the restriction of h on the ky axis, i.e. namely setting kx = 0, amounts to consider the particular case 
where species u and v are solely allowed to diffuse along the vertical direction. One can therefore equivalently set 
= 0 in Eq (|8]) and thus get: 


h{k,, ky) = det(J) - klifuD^y^ + ) + kp^y^D^y^ , (9) 

By solving equation (|9]) for ky^ one obtains two positive solutions, 0 < k- < if and only if the following conditions 
are met: 

> 0 

1 det(J) > 0. 


Let us observe that from the first relation of Eq. (Unil and the condition tr(J) < 0 implies for the 

instability to set in and the patterns to develop, the inhibitor should diffuse faster than the activator in the y direction. 

The symmetric limiting case is recovered when species u and v are allowed to diffuse only along the horizontal 
direction, which in turn amounts to restrict h to the kx axis. The analysis can be hence handled by setting = 
= 0 in Eq dHj) and proceeding in analogy with above. One can straightforwardly obtain the following necessary 
and sufficient conditions for the existence of Turing patterns: 

+ gM"^ > 0 , . 

+ gvD^D^ - det(J) > 0. 

Once again, from the first relation of Eq. m and the condition tr( J) < 0, one can immediately conclude that patterns 
are possible only if , namely if the inhibitor diffuses faster than the activator along the x direction. 

These conclusions are clearly not surprising, as they constitute an obvious adaptation of the standard Turing 
framework to the present context, in the trivial limit where one of the diffusion direction is alternatively silenced. 
Starting from this observation, it is however interesting to speculate on the possibility of turning unstable complex 
mixed modes (kx^ ky), via a symmetry breaking process of the Turing type, when the simplified pathways to pattern 
formation explored above are instead precluded. 

To this end, we go back to function h(kx^ky) and study its sign when moving on (kx^ky), along specific directions. 
More concretely, we set kx = ^ky, and vary the free parameter 7 to span the reference plane. Turing patterns can 
then develop only if h{^ky, ky) < 0 , where: 

h{^ky, ky) = det(J) - + g^D^y))] 

=; B^k^y-B2kl + Bs, ( 12 ) 

and the last expression defines the coefficients Bi , B 2 and Bs . It can be readily realized that Bi and Bs are positively 
definite, while B 2 can assume both positive and negative values. In the following, we shall impose the simultaneous 
violation of conditions m and m, via crossed negation of the corresponding inequalities, and look for possible 
values of the control parameter 7 that make the system unstable. 
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A. Conditions (llQD i and (llip i are not satisfied 


Let us thus assume 


f < 0 

< 0, 

while the remaining two conditions (p!Q|) ii and ([TT]) ii do hold. 

One can trivially realize that in this case B 2 is negative, hence h{'yky, ky) = Biky + \B 2 \ky B^ > 0 for all kx = "yky 
and ky. No instability can thus develop which seeds the emergence of self-organized Turing patterns. 


B. Conditions (ITOll i and ([TT]) ii are not satisfied 


We now assume 


ifuDir^+gM^^<0 

det(J) < 0, 

while the remaining two relations (fTQ|) ii and (pT]) i are verified. 

Solving for the limiting condition h{jky, ky) = 0 one gets a closed expression for ky. By imposing ky to be positive 
yields ^2 > 0 and B 2 — ABiB^, > 0 . 

A straightforward computation gives: 

^2 > 0 if 7 ^ > gi, 

where 

=-fJl-’ 

fuD^v^ + 

where use has been made of Eq. (dlDi. 

A somehow lengthy computation allows us to write: 

Bl - AB 1 B 3 = Ai7^ + ^ 27 " + ^3 , (13) 


where: 

El-4det(J)D7ll7 (14) 

2(/„7") + g,Dl^^){UDiy^ + g.D^y^) - 4det(J)(!?(*')(15) 

r 2 -4det(J)E)yD7 (16) 

ri = (/„D7+fl„il7)2 and r 2 = (/„D7+fl„D7)2, (17) 

Under the above assumptions Fi < Adet{J)D^^Di^\ which implies Ai < 0. Similarly, as r 2 > 4 det( 

^3 > 0 . On the other hand, A 2 < 0, this latter quantity resulting from the sum of two negative terms. Hence, 

B 2 — 4 H 1 H 3 > 0 if 0 < 7 ^ < g 2 , where q 2 = ^ q 

We can easily show that qi > g 2 , which in turn implies that B 2 and — 4 H 1 H 3 cannot be at the same time positive, 
as it should happen for the instability to develop. We can hence conclude that Turing patterns cannot develop in this 
case either. 


= 

A2 = 
^3 = 
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C. Conditions (llOp ii and (llip i are not satisfied 


Let us thus assume 

I det(J) < 0 

< 0 , 

while the remaining two condition([T0|)i and (pT]) ii are verified. 

Once again requiring h{jky,ky) < 0 , necessarily imply ^2 > 0 and > 0 . The former condition is 

satisfied whenever: 

7^ e (0,gi), 

for qi = —{fuD^v^ + g-uD^^) / {fuD^v^ + g^D^^) > 0 . The latter condition — 4i?ii?3 > 0 yields 

7 " (Ti - 4det(J)Dl")Z)(")) + 7 '[ 2 (/„Z)(-) + + g^D^y^) - 4det(J)(Z)(*^)D(*) + 

+ Ta - 4 det{J)Dly^Dly'> ;= jUi + j^A2 + A 3 > 0. 

Here, H i > 0 while H2 < 0 and H3 < 0 . Hence, the previous inequality is satisfied for any 7^ > q2 for q2 = 
{—A2 + ^/H2'^^T4iH3)/(2Hi) > 0 . However, one can prove that qi < g2, which implies that B2 and — dHiHs 
cannot be simultaneously positive. The conclusion is therefore that h{'yky, ky) > 0 , and Turing patterns cannot take 
place. 


D. Conditions (ITOll ii and ([TT]) ii are not satisfied 


Let us thus assume 

f {fuD^^^ + gM^^? - det(J) < 0 

det(J) < 0, 

while the remaining two assumptions (pT|) i and (EJi do hold. 

Under the present working hypothesis, the coefficient Hi, B2 and H3 are positive. Thus h{'yky, ky) can take negative 
values, if and only if H| — ABiB^ > 0. As previously remarked, we can rewrite 

Bl - 4B1B3 = A174 + A27" + As , 

where for i = 1, 2, 3 are defined as in (p! 4 |) . One can show that Hi and H3 are negative while H2 can take both signs. 
To satisfy the requirement H| — 4H1H3 > 0 the conditions H2 > 0 and — 4H1H3 > 0 should be simultaneously 
met. 

Let us rewrite H2 as follows 

A2 = 2711 - 4det(J) + Di^^Diy^) , 

where T^ have been defined in Eq. (ED. Straightforward manipulations allow us to write: 


A2 = 27l\7fI-4det(J) 

v ^ 


{x) 


(y) 


< — [ r2 


D 


{x) 


D 


{y) 


Ti- 


D 


{y )' 


D 


{x) 


]j(x) jj{x) 


p^\ 

^ 7 ; 



< 0. 


Since H2 is bound to be negative, the condition for Turing instability h{jky, ky) < 0 cannot be satisfied. 

Summing up we have demonstrated that patterns can eventually develop only if the system can undergo a symmetry 
breaking instability of the Turing type, in its restricted configuration where the diffusion is solely allowed along one 
spatial direction, either x or y. The result is summarized in Figure [U where different types of instabilities are 
schematically depicted. 

Interestingly, the instability can set in also if the inhibitor diffuses slower that the activator along one selected 
direction, provided the opposite holds for the transport along the orthogonal direction. In this respect, accounting for 
anisotropic diffusion enables one to partially relax the stringent conditions that underly the formation of the Turing 
motifs. In the next section, we will built on this observation and provide a numerical demonstration of the investigated 
phenomenon. 
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case b 




case c 


case d 


FIG. 1: Possible types of instabilities. Case a: h{kx,ky) > 0 for all > 0 and ky > 0. The system cannot turn unstable. 
Case b: h restricted to the kx axis takes negative values: a bounded contiguous domain in kx > 0 and ky>0 exist, for which 
h(kx^ky) > 0. Case c: h restricted to the ky axis takes negative values. Again a portion of the reference plan, adjacent to 
the domain of instability in kx = 0, can be found where h{kx, ky) > 0: Case d: the system is unstable along both kx = 0 and 
ky = 0 directions. The instability also interests non trivial modes with both kx ^ 0 and ky 7^= 0. 


IV. NUMERICAL ANALYSIS 


The aim of this section is to discuss a numerical implementation of the theory presented above. In particular, 
we will show that complex patterns can emerge for a system of two species in mutual interaction and undergoing 
anisotropic diffusion, also if the conventional Turing request of having inhibitors faster than activators is relaxed, 
along one of the two orthogonal directions of movements. To perform the analysis we operate in the framework of the 
so called Mimura-Murray model [H. The quantities u and v can be associated to prey and predator densities, which 
interact via the non-linear functions: 


f{u, v) = ((a -\-hu — v?)lc — v) u and g{u^ 'T’) = — (1 + dv)) v ; (18) 

the model possesses 6 equilibria, whose stability and positivity depend on the value of the chosen parameters. We 
here focus on the fixed point {u, v) 


u = 1 -\- 


bd — 2d — c-\- y/~K 


2 d 


and V = 


bd — 2d — c -\- 




where A = {bd — 2d — c)^ + 4d^(a + 6 — 1), (19) 


and assume a = 35, 5 = 16, c = 9 and d = 0.4 which in turn implies {u^v) = (5,10). Moreover, the Jacobian entries 
evaluated at the fixed point reads fu = 3.33, fv = —5, gu = 10 and gy = —4. Hence, det(J) > 0 and tr(J) < 0: 
the fixed point is a stable equilibrium. We also remark that u acts as the activator and v stands for the inhibitor 
species, as /^ > 0 and gy < 0. Under specific conditions, the fixed point can be destabilized by an external, non 
homogeneous, perturbation, paving the way to the subsequent generation of Turing patterns, in the non linear regime 
of the evolution. In Fig.[2]we report a gallery of representative patterns that can be obtained under distinct conditions. 

To generate the asymptotic patterns displayed in panel (a) of Fig[2l parameters are set so that both relations (pT|) 
and dm) are satisfied, > D^\c and > Du^Vc, where Vy ^ 16. Inhibitor diffuses faster than activators 

in both X and y directions, although with different diffusion constants. The dispersion relation (see Fig[2fb)) can 
be assimilated to that sketched in Fig. [Hd), and the corresponding patterns share marked similarities with those 
obtained in the conventional case of isotropic transport. 

In panel (c) of Fig. O conditions (pT|) hold, while ([TT]) do not, while < D^\c, where Vy ^ 16. 

The dispersion relation. Fig. [2fd), is also depicted and shown to resemble that displayed in Fig. Hl^c). The patterns 
which follow this unusual choice of the diffusion constants, compare nicely with those emerging under the standard 
paradigm, this is because D^y^/D^u^ is smaller but close to Vy. 

Finally, in panel (e,f) of FigO the activator is assigned a diffusion coefficient Du^ is larger than D^\ the homologous 

(x) (x) 

constant associated to the inhibitor species, and still Dy ^ > Du Vy. The dispersion relation falls in the category 
exemplified in Fig. [U^c), and the corresponding patterns are found to organize in regular stripes, which run almost 
parallel to the direction where the instability is present. 
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FIG. 2: Asymptotic activator distribution: the concentration u{x^ y, t) is displayed for sufficiently large t. Panel (a): — 0.01, 

= 0.16, = 0.02, = 0.32. Panel (b): = 0.01, = 0.155, = 0.02, = 0.32. Panel (c): 

— 0.012, — 0.01, — 0.02, — 0.32. The other parameters are set to a = 35; 5 = 16; c = 9; d = 0.4. The hlled 

black squares identify the position of the maximum of the dispersion relation. 


V. CONCLUSIONS 

In this paper we elaborated on the impact of anisotropic diffusion for the emergence of Turing patterns in reaction- 
diffusion systems. We have in particular focused on systems of two interacting species confined in a rectangular, 
continuum domain, endowed with periodic boundary conditions. With reference to this paradigmatic case study, we 
have shown that a symmetry breaking instability of the Turing type can occur only if patterns do exist when diffusion 
is impeded along one of the two accessible directions. In other words, patterns which resemble those obtained in 
the conventional setting of isotropic diffusion emerge, also when the standard Turing condition is violated along one 
specific direction. Interestingly, the instability can also occur if the activator diffuses faster than the inhibitor, along 
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the direction of spatial relocation for which the usual Turing conditions are not met. 
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